/*
  File:    extend.cpp
  Purpose: These routines extend and contract matrices by several layers. 
           Useful for routines like WENO, that may require more points than are
           available.

  Contact: Paul Macklin
           pmacklin@math.uci.edu
		  http://math.uci.edu/~pmacklin
*/

#include <cstdlib>
#include <cstdio>
#include <cmath>

using namespace std;

#include "../Matrix.h"

Matrix extend(Matrix &mIn, int nLayers)
{
 Matrix mOut(mIn.TellRows()+2*nLayers, mIn.TellCols()+2*nLayers);
 int i=0,j=0;

 for(j=1+nLayers; j <= mIn.TellCols()+nLayers; j++)
 {
  double dSlopeXBottom = 
         (*mIn(2-1,j-1-nLayers)-*mIn(1-1,j-1-nLayers)+0.0); 
  double dSlopeXTop = 
         (*mIn(mIn.TellRows()-1,j-1-nLayers) - 
          *mIn(mIn.TellRows()-1-1,j-1-nLayers)+0.0); 
  for(i=1; i <= nLayers; i++)
  {      
   *mOut(i-1,j-1) = *mIn(1-1,j-1-nLayers) - 
          (nLayers-i+1)*dSlopeXBottom + 0.0;
  }
  for(i=nLayers+1; i <= nLayers+mIn.TellRows() ; i++)
  {
   *mOut(i-1,j-1) = *mIn(i-1-nLayers,j-1-nLayers);
  }
  for(i=nLayers+mIn.TellRows()+1; i <= nLayers+mIn.TellRows()+nLayers; i++)
  {
   *mOut(i-1,j-1) = *mIn(mIn.TellRows()-1,j-1-nLayers) + 
          (i - (nLayers+mIn.TellRows()))*dSlopeXTop;
  }
 }

 for(i=1; i<= mIn.TellRows()+2*nLayers ; i++)
 {
  double dSlopeXLeft=(*mOut(i-1,nLayers+2-1)-
                      *mOut(i-1,nLayers+1-1));
  double dSlopeXRight=(*mOut(i-1,nLayers-1+mIn.TellCols())-
                       *mOut(i-1,nLayers-1+mIn.TellCols()-1));
  for(j=1 ; j <= nLayers ; j++)
  {
   *mOut(i-1,j-1) = *mOut(i-1,nLayers-1+1)-
                      (nLayers-j+1)*dSlopeXLeft;
  }
  for(j=nLayers+mIn.TellCols()+1; j <= nLayers+mIn.TellCols()+nLayers ; j++)
  {
   *mOut(i-1,j-1) = *mOut(i-1,nLayers-1+mIn.TellCols())+
                      (j-(nLayers+mIn.TellCols()))*dSlopeXRight;
  }
 }
 return mOut;
}

Matrix contract(Matrix &mIn, int nLayers)
{
 Matrix mOut = zeros(mIn.TellRows()-2*nLayers, mIn.TellCols()-2*nLayers);
 int i=0, j=0;
 for(i=0; i < mOut.TellRows() ; i++)
 {for(j=0; j < mOut.TellCols() ; j++)
  {
   *mOut(i,j) = *mIn(i+nLayers,j+nLayers);
  }
 }
 return mOut; 
}
